#ifndef __MY_MATRIX_H__
#define __MY_MATRIX_H__

#include "defines.h"
#include <iostream>
#include <math.h>

// https://blog.csdn.net/darkrabbit/article/details/80025935
// https://blog.csdn.net/a133900029/article/details/80445198
template <class T>
class Matrix2D {
private:
	using val_t = T;
	static constexpr val_t EPS = val_t(1e-6);

	int row_, col_;
	val_t *data_;

	static bool equals(val_t a, val_t b) { return fabs(a - b) <= EPS; }
	static bool zero(val_t a) { return fabs(a) <= EPS; }

public:
	Matrix2D(int row = 0, int col = 0, val_t *data = nullptr) : data_(nullptr) {
		if (row > 0 && col > 0) {
			_create(row, col);
			if (data) memcpy(data_, data, row * col * sizeof(val_t));
		}
	}

	Matrix2D(const Matrix2D& m) {
		if (m.row_ > 0 && m.col_ > 0) {
			_create(m.row_, m.col_);
			if (m.data_) memcpy(data_, m.data_, m.size() * sizeof(val_t));
		}
	}

	~Matrix2D() { _empty(); }

	int row() const { return row_; }
	int col() const { return col_; }
	val_t * data() const { return data_; }

private:
	void _empty() { SAFE_DELETE_ARR(data_); }
	bool _sameSize(const Matrix2D& m) const { return row_ == m.row_ && col_ == m.col_; }
	void _create(int row, int col) {
		int size = row * col;
		data_ = new val_t[size];
		memset(data_, 0, size * sizeof(val_t));
		row_ = row; col_ = col;
	}

	// 删除第r行，第c列后的矩阵
	Matrix2D _subMatrix(int r, int c) const {
		Matrix2D _tmp(row_ - 1, col_ - 1);
		int pos = 0;
		for (int y = 0, offset = 0; y < row_; y++, offset += col_) {
			if (y == r) continue;
			for (int x = 0; x < col_; x++) {
				if (x == c) continue;
				_tmp.data_[pos++] = data_[offset + x];
			}
		}
		return _tmp;
	}

	// https://www.shuxuele.com/algebra/matrix-determinant.html
	val_t _det() const {
		if (row_ == 1) return data_[0];
		if (row_ == 2) return data_[0] * data_[3] - data_[1] * data_[2];

		val_t res = 0;
		int flag = 1;
		for (int i = 0; i < col_; i++) {
			res += flag * data_[i] * _subMatrix(0, i)._det();
			flag *= -1;
		}
		return res;
	}

	val_t _tr() const {
		val_t sum = 0;
		for (int m = 0, offset = 0; m < row_; m++, offset += col_) {
			sum += data_[offset + m];
		}
		return sum;
	}

	// https://blog.csdn.net/beauthy/article/details/117422074
	Matrix2D _adj() const {
		Matrix2D _tmp(row_, col_);
		for (int y = 0, offset = 0; y < _tmp.row_; y++, offset += _tmp.col_) {
			for (int x = 0; x < _tmp.col_; x++) {
				auto m = std::move(_subMatrix(y, x));
				auto detm = m._det();
				int sign = is_odd(x + y) ? -1 : 1;
				_tmp.data_[offset + x] = sign * detm;
			}
		}
		return _tmp.t();
	}

	// https://jingyan.baidu.com/article/925f8cb8a74919c0dde056e7.html
	Matrix2D _inv() const {
		val_t res = _det();
		if (zero(res)) return Matrix2D();
		float ratio = 1.0f / res;
		return ratio * _adj();
	}

public:
	bool isEmpty() const { return data_ != nullptr; }
	int size() const { return row_ * col_; }

	// 转置 transpose
	Matrix2D t() const {
		Matrix2D _tmp(col_, row_);
		for (int y = 0, offset = 0; y < _tmp.row_; y++, offset += _tmp.col_) {
			for (int x = 0; x < _tmp.col_; x++) {
				_tmp.data_[offset + x] = data_[x * col_ + y];
			}
		}
		return _tmp;
	}

	// 伴随矩阵 adjugate
	Matrix2D adj() const {
		if (row_ != col_) return Matrix2D();
		return _adj();
	}

	// 逆 Inverse  A*/|A|
	Matrix2D inv() const {
		if (row_ != col_) return Matrix2D();
		return _inv();
	}

	// 迹(Trace) 方阵A的对角线之和称为迹，记作tr(A)
	bool tr(val_t& /*out*/ res) const {
		if (row_ != col_) return false;
		res = _tr();
		return true;
	}

	// 行列式 |m|. 只有方阵才能定义行列式
	bool det(val_t& /**/res) const {
		if (row_ != col_) return false;
		res = _det();
		return true;
	}

	friend std::ostream& operator << (std::ostream &os, const Matrix2D& m) {
		for (int y = 0, offset = 0; y < m.row_; y++, offset += m.col_) {
			for (int x = 0; x < m.col_; x++) {
				os << m.data_[offset + x] << "\t";
			}
			os << std::endl;
		}
		return os;
	}

	val_t & operator [] (int i) {
		return data_[i];
	}

	Matrix2D& operator = (const Matrix2D& m) {
		if (!_sameSize(m)) {
			_empty();
			_create(m.row_, m.col_);
		}
		if (m.data_) memcpy(data_, m.data_, row_ * col_ * sizeof(val_t));
		return *this;
	}

	Matrix2D operator + (const Matrix2D& m) const {
		if (!_sameSize(m)) return Matrix2D(); // empty matrix
		Matrix2D _tmp(row_, col_, data_);
		for (int i = 0, N = _tmp.size(); i < N; i++) {
			_tmp.data_[i] += m.data_[i];
		}
		return _tmp;
	}

	Matrix2D operator - (const Matrix2D& m) const {
		if (!_sameSize(m)) return Matrix2D(); // empty matrix
		Matrix2D _tmp(row_, col_, data_);
		for (int i = 0, N = _tmp.size(); i < N; i++) {
			_tmp.data_[i] -= m.data_[i];
		}
		return _tmp;
	}

	Matrix2D operator * (const Matrix2D& m) const {
		if (col_ != m.row_) return Matrix2D();
		Matrix2D _tmp(row_, m.col_);
		for (int y = 0, offset = 0; y < _tmp.row_; y++, offset += _tmp.col_) {
			for (int x = 0; x < _tmp.col_; x++) {
				val_t sum = 0;
				for (int ls = y * col_, le = ls + col_, rs = x; ls < le; ls++, rs += m.col_) {
					sum += data_[ls] * m.data_[rs];
				}
				_tmp.data_[offset + x] = sum;
			}
		}
		return _tmp;
	}

	Matrix2D operator * (val_t a) const {
		Matrix2D _tmp(row_, col_);
		int offset = 0;
		for (int y = 0; y < _tmp.row_; y++) {
			for (int x = 0; x < _tmp.col_; x++) {
				int pos = offset + x;
				_tmp.data_[pos] = data_[pos] * a;
			}
			offset += _tmp.col_;
		}
		return _tmp;
	}

	bool operator == (const Matrix2D& m) const { return !operator!=(m); }
	bool operator != (const Matrix2D& m) const {
		if (!_sameSize(m)) return true;
		for (int i = 0, N = size(); i < N; i++) {
			if (!equals(data_[i], m.data_[i])) return true;
		}
		return false;
	}

	friend Matrix2D operator * (val_t a, const Matrix2D& m) {
		return m.operator*(a);
	}
};

typedef Matrix2D<float> Matrix2Df;

#endif // __MY_MATRIX_H__